Introduction to NeXL¶

Nicholas W. M. Ritchie 26-Apr-2023 (nicholas.ritchie@nist.gov)

What is NeXL?¶

NeXL is a collection of packages (libraries) for the Julia language that implement common X-ray microanalysis algorithms, physical data and utilities.

The core NeXL packages are

  • NeXLCore - Defines core concepts like elements, X-rays, atomic shells, k-ratios and fundamental electron and X-ray algorithms.
  • NeXLMatrixCorrection - Implements algorithms for performing matrix correction for bulk samples
  • NeXLSpectrum - Implements data types and algorithms to handle spectra and hyper-spectra

NeXL builds on many packages from the Julia infrastructure including

  • Gadfly - A "grammar of graphics" plotting package
  • DataFrames - "Pandas-like" tools for working with tabular data

Why Julia?¶

Julia is a high-performance scripting language designed for numerical data analysis. The base Julia libraries provide basic scalar, vector and matrix numerical analysis tools. Julia differs from most scriptable languages in that all code is compiled before execution. This leads to a "slow then fast" user experience. The first time a method is used, it and all the methods it depends upon must be compiled. This can lead to a significant delay. However, once compiled, the code is fast - often as fast as C or Fortran code. You will notice this "slow then fast" nature as you work with this notebook. It is particularly evident when creating plots.

NeXLCore¶

We will start by exploring NeXLCore since this package provides core functionality for the other NeXL packages.

In [ ]:
using DrWatson
@quickactivate "IUMAS"

using NeXLCore
  Activating project at `c:\Users\nritchie\Desktop\IUMAS`

Elements¶

We will first construct a datum representing an element and then manipulate this object to extract elemental data.

Julia is not an object-oriented language. Instead it uses what is called "multiple dispatch" to determine which method implementation to call dependent on the arguments handed to the function.

In [ ]:
el = n"Tc"  # Here we are using a special syntactic sugar to convert the string "Tc" into an Element struct
Out[ ]:
Technetium (Tc), number 43:
categorytransition metal
atomic mass98.0 u
density11.0 g/cm³
molar heat24.27 J/mol⋅K
melting point2430.0 K
boiling point4538.0 K
phaseSolid
shells[2, 8, 18, 13, 2]
electron configuration1s² 2s² 2p⁶ 3s² 3p⁶ 4s² 3d¹⁰ 4p⁶ 5s² 4d⁵
appearanceshiny gray metal
summaryTechnetium (/tɛkˈniːʃiəm/) is a chemical element with symbol Tc and atomic number 43. It is the element with the lowest atomic number in the periodic table that has no stable isotopes:every form of it is radioactive. Nearly all technetium is produced synthetically, and only minute amounts are found in nature.
discovered byEmilio Segrè
sourcehttps://en.wikipedia.org/wiki/Technetium

Let's get some properties of this elemental datum.

In [ ]:
z(el), a(el), symbol(el), name(el)
Out[ ]:
(43, 98.0, "Tc", "Technetium")

We can combine elements into materials. The mat"..." macro converts expressions into Material structs.

In [ ]:
mat = mat"Ca(PO4)3OH"
mat, typeof(mat)
Out[ ]:
(Ca(PO4)3OH[P=0.2717,Ca=0.1172,O=0.6082,H=0.0029], Material{Float64, Float64})

An alternative syntax allows you to specify additional properties of the Material.

In [ ]:
mat = parse(Material, "Ca(PO4)3OH", name = "Apatite", density = 2.2)
Out[ ]:
Apatite[P=0.2717,Ca=0.1172,O=0.6082,H=0.0029,2.20 g/cm³]

In addition to being able to parse chemical formulae, NeXLCore can also parse expressions like the next one. In this case, the multiplier represents the mass fraction of the element. The right hand side can also be a chemical formula as though the material was made from various mass fractions of consituent materials.

In [ ]:
adm = parse(Material, "0.3399*O+0.0664*Al+0.0405*Si+0.0683*Ca+0.0713*Ti+0.1055*Zn+0.3037*Ge", name="ADM-6005a")
Out[ ]:
ADM-6005a[Ti=0.0713,Al=0.0664,Ge=0.3037,Ca=0.0683,Zn=0.1055,O=0.3399,Si=0.0405]

In addition to elements and materials, NeXLCore provides mechanisms to work with characteristic X-rays, atomic shells and other atomic physics data. The n"..." macro can be used to create Element, CharXRay, AtomicSubShell and SubShell structs.

In [ ]:
cxr = n"Tc K-L3"
cxr, typeof(cxr)
Out[ ]:
(Tc K-L3, CharXRay)
In [ ]:
ass = n"Tc K"
ass, typeof(ass)
Out[ ]:
(Tc K, AtomicSubShell)

Let's get some properties of the CharXRay datum. You'll notice that the energy is in electron-volts as is the case for all energies within NeXL.

In [ ]:
energy(cxr), inner(cxr), energy(inner(cxr)), outer(cxr), energy(outer(cxr)), typeof(inner(cxr))
Out[ ]:
(18367.0, Tc K, 21050.0, Tc L3, 2683.0, AtomicSubShell)

To determine every characteristic X-ray generated by a particular element, you can use the characteristic(...). To generate sub-sets of the transitions replace alltransitions with ktransitions, ltransitions or mtransitions.

In [ ]:
cxrs=characteristic(el, alltransitions)
cxrs, length(cxrs)
Out[ ]:
(Tc K-L3 + 9 others, Tc L3-M5 + 23 others & Tc M5-N3 + 9 others, 44)

Now we can use various functions to extract X-ray related physical data from elements and materials. The mac(...) function returns the mass absorption coeffficient in $cm^2/g$.

In [ ]:
mac(mat, n"Ca K-L3"), mac(n"Ca", n"Ca K-L3"), mac(mat, 3692.0), mac(n"Ca", 3692.0)
Out[ ]:
(257.51384578832244, 139.20100028283625, 257.51384578832244, 139.20100028283625)

This is an example of one of Julia's core design features - multiple dispatch. The function mac(...) is called on four distinct argument types.

  • mac(Material, CharXRay)
  • mac(Element, CharXRay)
  • mac(Material, Float64)
  • mac(Element, Float64)

This is similar but extends on the way object oriented code can use the object type to determine which function implementation to call.

NeXL uses the Gadfly library to plot spectra and other data items.

In [ ]:
using Gadfly
plot(e->mac(mat, e), 100.0, 1.0e4, Scale.y_log10)
Out[ ]:
x -1.50×10⁴ -1.00×10⁴ -5.00×10³ 0 5.00×10³ 1.00×10⁴ 1.50×10⁴ 2.00×10⁴ 2.50×10⁴ -1.00×10⁴ -9.00×10³ -8.00×10³ -7.00×10³ -6.00×10³ -5.00×10³ -4.00×10³ -3.00×10³ -2.00×10³ -1.00×10³ 0 1.00×10³ 2.00×10³ 3.00×10³ 4.00×10³ 5.00×10³ 6.00×10³ 7.00×10³ 8.00×10³ 9.00×10³ 1.00×10⁴ 1.10×10⁴ 1.20×10⁴ 1.30×10⁴ 1.40×10⁴ 1.50×10⁴ 1.60×10⁴ 1.70×10⁴ 1.80×10⁴ 1.90×10⁴ 2.00×10⁴ -1.0×10⁴ 0 1.0×10⁴ 2.0×10⁴ -1.000×10⁴ -9.500×10³ -9.000×10³ -8.500×10³ -8.000×10³ -7.500×10³ -7.000×10³ -6.500×10³ -6.000×10³ -5.500×10³ -5.000×10³ -4.500×10³ -4.000×10³ -3.500×10³ -3.000×10³ -2.500×10³ -2.000×10³ -1.500×10³ -1.000×10³ -5.000×10² 0 5.000×10² 1.000×10³ 1.500×10³ 2.000×10³ 2.500×10³ 3.000×10³ 3.500×10³ 4.000×10³ 4.500×10³ 5.000×10³ 5.500×10³ 6.000×10³ 6.500×10³ 7.000×10³ 7.500×10³ 8.000×10³ 8.500×10³ 9.000×10³ 9.500×10³ 1.000×10⁴ 1.050×10⁴ 1.100×10⁴ 1.150×10⁴ 1.200×10⁴ 1.250×10⁴ 1.300×10⁴ 1.350×10⁴ 1.400×10⁴ 1.450×10⁴ 1.500×10⁴ 1.550×10⁴ 1.600×10⁴ 1.650×10⁴ 1.700×10⁴ 1.750×10⁴ 1.800×10⁴ 1.850×10⁴ 1.900×10⁴ 1.950×10⁴ 2.000×10⁴ h,j,k,l,arrows,drag to pan i,o,+,-,scroll,shift-drag to zoom r,dbl-click to reset c for coordinates ? for help ? 10-6 10-5 10-4 10-3 10-2 10-1 100 101 102 103 104 105 106 107 108 109 1010 1011 10-5.0 10-4.5 10-4.0 10-3.5 10-3.0 10-2.5 10-2.0 10-1.5 10-1.0 10-0.5 100.0 100.5 101.0 101.5 102.0 102.5 103.0 103.5 104.0 104.5 105.0 105.5 106.0 106.5 107.0 107.5 108.0 108.5 109.0 109.5 1010.0 10-5 100 105 1010 10-5.0 10-4.8 10-4.6 10-4.4 10-4.2 10-4.0 10-3.8 10-3.6 10-3.4 10-3.2 10-3.0 10-2.8 10-2.6 10-2.4 10-2.2 10-2.0 10-1.8 10-1.6 10-1.4 10-1.2 10-1.0 10-0.8 10-0.6 10-0.4 10-0.2 100.0 100.2 100.4 100.6 100.8 101.0 101.2 101.4 101.6 101.8 102.0 102.2 102.4 102.6 102.8 103.0 103.2 103.4 103.6 103.8 104.0 104.2 104.4 104.6 104.8 105.0 105.2 105.4 105.6 105.8 106.0 106.2 106.4 106.6 106.8 107.0 107.2 107.4 107.6 107.8 108.0 108.2 108.4 108.6 108.8 109.0 109.2 109.4 109.6 109.8 1010.0 f(x)
In [ ]:
plot(
    layer(e->mac(mat, e), 100.0, 1.0e4, Theme(default_color="Red")),
    layer(e->mac(n"Ca", e), 100.0, 1.0e4, Theme(default_color="Green")),
    layer(e->mac(n"P", e), 100.0, 1.0e4, Theme(default_color="Blue")),
    Scale.y_log10
)
Out[ ]:
xmin -1.50×10⁴ -1.00×10⁴ -5.00×10³ 0 5.00×10³ 1.00×10⁴ 1.50×10⁴ 2.00×10⁴ 2.50×10⁴ -1.00×10⁴ -9.00×10³ -8.00×10³ -7.00×10³ -6.00×10³ -5.00×10³ -4.00×10³ -3.00×10³ -2.00×10³ -1.00×10³ 0 1.00×10³ 2.00×10³ 3.00×10³ 4.00×10³ 5.00×10³ 6.00×10³ 7.00×10³ 8.00×10³ 9.00×10³ 1.00×10⁴ 1.10×10⁴ 1.20×10⁴ 1.30×10⁴ 1.40×10⁴ 1.50×10⁴ 1.60×10⁴ 1.70×10⁴ 1.80×10⁴ 1.90×10⁴ 2.00×10⁴ -1.0×10⁴ 0 1.0×10⁴ 2.0×10⁴ -1.000×10⁴ -9.500×10³ -9.000×10³ -8.500×10³ -8.000×10³ -7.500×10³ -7.000×10³ -6.500×10³ -6.000×10³ -5.500×10³ -5.000×10³ -4.500×10³ -4.000×10³ -3.500×10³ -3.000×10³ -2.500×10³ -2.000×10³ -1.500×10³ -1.000×10³ -5.000×10² 0 5.000×10² 1.000×10³ 1.500×10³ 2.000×10³ 2.500×10³ 3.000×10³ 3.500×10³ 4.000×10³ 4.500×10³ 5.000×10³ 5.500×10³ 6.000×10³ 6.500×10³ 7.000×10³ 7.500×10³ 8.000×10³ 8.500×10³ 9.000×10³ 9.500×10³ 1.000×10⁴ 1.050×10⁴ 1.100×10⁴ 1.150×10⁴ 1.200×10⁴ 1.250×10⁴ 1.300×10⁴ 1.350×10⁴ 1.400×10⁴ 1.450×10⁴ 1.500×10⁴ 1.550×10⁴ 1.600×10⁴ 1.650×10⁴ 1.700×10⁴ 1.750×10⁴ 1.800×10⁴ 1.850×10⁴ 1.900×10⁴ 1.950×10⁴ 2.000×10⁴ h,j,k,l,arrows,drag to pan i,o,+,-,scroll,shift-drag to zoom r,dbl-click to reset c for coordinates ? for help ? 10-7 10-6 10-5 10-4 10-3 10-2 10-1 100 101 102 103 104 105 106 107 108 109 1010 1011 1012 1013 10-6.0 10-5.5 10-5.0 10-4.5 10-4.0 10-3.5 10-3.0 10-2.5 10-2.0 10-1.5 10-1.0 10-0.5 100.0 100.5 101.0 101.5 102.0 102.5 103.0 103.5 104.0 104.5 105.0 105.5 106.0 106.5 107.0 107.5 108.0 108.5 109.0 109.5 1010.0 1010.5 1011.0 1011.5 1012.0 10-10 100 1010 1020 10-6.0 10-5.8 10-5.6 10-5.4 10-5.2 10-5.0 10-4.8 10-4.6 10-4.4 10-4.2 10-4.0 10-3.8 10-3.6 10-3.4 10-3.2 10-3.0 10-2.8 10-2.6 10-2.4 10-2.2 10-2.0 10-1.8 10-1.6 10-1.4 10-1.2 10-1.0 10-0.8 10-0.6 10-0.4 10-0.2 100.0 100.2 100.4 100.6 100.8 101.0 101.2 101.4 101.6 101.8 102.0 102.2 102.4 102.6 102.8 103.0 103.2 103.4 103.6 103.8 104.0 104.2 104.4 104.6 104.8 105.0 105.2 105.4 105.6 105.8 106.0 106.2 106.4 106.6 106.8 107.0 107.2 107.4 107.6 107.8 108.0 108.2 108.4 108.6 108.8 109.0 109.2 109.4 109.6 109.8 1010.0 1010.2 1010.4 1010.6 1010.8 1011.0 1011.2 1011.4 1011.6 1011.8 1012.0 y

There is much more functionality in NeXLCore, but that is enought to get us started. Next let's take a look at NeXLSpectrum.

NeXLSpectrum¶

NeXLSpectrum defines methods for reading, writing, manipulating, plotting and fitting X-ray spectra and X-ray spectrum images ("hyperspectra").

In [ ]:
using NeXLSpectrum
[ Info: Loading Gadfly support into NeXLSpectrum.

First, we will read a spectrum from disk using the loadspectrum(...) method. It take a filename and is able to sniff the file type for many common spectrum file formats.

In [ ]:
sp = loadspectrum(joinpath(datadir(), "exp_raw", "ADM6005a spectra","ADM-6005a_1.msa"))
Out[ ]:
  *                                                                                                  128860.0
  *                                                                                                 
  *                                                                                                 
  *                                                                                                 
  *                                                                                                 
  *                                                                                                 
  *     *                                                                                           
  *     *                                                                                           
  *     *                                                                                           
  * ** **                                                                                           
  * ** *****        **  **                                                                          
**************************************************************************************************** 20.0 keV
Spectrum{Float64}[ADM-6005a_1, -484.20818 + 5.01716⋅ch eV, 4096 ch, 20.0 keV, unknown composition, 6.81e6 counts]
In [ ]:
set_default_plot_size(10inch, 4inch)
spec_files = [ joinpath(datadir(), "exp_raw", "ADM6005a spectra","ADM-6005a_$i.msa") for i in 1:15 ]
specs=loadspectrum.(spec_files)
plot(specs..., klms=[n"O", n"Al", n"Si", n"Ca", n"Ti", n"Zn", n"Ge"])
Out[ ]:
Energy (eV) -25000 -20000 -15000 -10000 -5000 0 5000 10000 15000 20000 25000 30000 35000 40000 45000 -20000 -18000 -16000 -14000 -12000 -10000 -8000 -6000 -4000 -2000 0 2000 4000 6000 8000 10000 12000 14000 16000 18000 20000 22000 24000 26000 28000 30000 32000 34000 36000 38000 40000 -20000 0 20000 40000 -20000 -19000 -18000 -17000 -16000 -15000 -14000 -13000 -12000 -11000 -10000 -9000 -8000 -7000 -6000 -5000 -4000 -3000 -2000 -1000 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 11000 12000 13000 14000 15000 16000 17000 18000 19000 20000 21000 22000 23000 24000 25000 26000 27000 28000 29000 30000 31000 32000 33000 34000 35000 36000 37000 38000 39000 40000 Si Si Ti Ti Ti Ti Ge Ge Al Al Ca Ca Ca Zn Zn Zn O O Ti Ti Ti Ca Ca Ca Ca Ge Ge Ge Zn Zn h,j,k,l,arrows,drag to pan i,o,+,-,scroll,shift-drag to zoom r,dbl-click to reset c for coordinates ? for help ? ADM-6005a_1 ADM-6005a_2 ADM-6005a_3 ADM-6005a_4 ADM-6005a_5 ADM-6005a_6 ADM-6005a_7 ADM-6005a_8 ADM-6005a_9 ADM-6005a_10 ADM-6005a_11 ADM-6005a_12 ADM-6005a_13 ADM-6005a_14 ADM-6005a_15 Spectra -100000 -80000 -60000 -40000 -20000 0 20000 40000 60000 80000 100000 120000 140000 160000 180000 -65000 -60000 -55000 -50000 -45000 -40000 -35000 -30000 -25000 -20000 -15000 -10000 -5000 0 5000 10000 15000 20000 25000 30000 35000 40000 45000 50000 55000 60000 65000 70000 75000 80000 85000 90000 95000 100000 105000 110000 115000 120000 125000 -100000 0 100000 200000 -64000 -62000 -60000 -58000 -56000 -54000 -52000 -50000 -48000 -46000 -44000 -42000 -40000 -38000 -36000 -34000 -32000 -30000 -28000 -26000 -24000 -22000 -20000 -18000 -16000 -14000 -12000 -10000 -8000 -6000 -4000 -2000 0 2000 4000 6000 8000 10000 12000 14000 16000 18000 20000 22000 24000 26000 28000 30000 32000 34000 36000 38000 40000 42000 44000 46000 48000 50000 52000 54000 56000 58000 60000 62000 64000 66000 68000 70000 72000 74000 76000 78000 80000 82000 84000 86000 88000 90000 92000 94000 96000 98000 100000 102000 104000 106000 108000 110000 112000 114000 116000 118000 120000 122000 124000 126000 Counts
In [ ]:
plot(specs..., klms=[n"O", n"Al", n"Si", n"Ca", n"Ti", n"Zn", n"Ge"], xmax=6.0e3)
Out[ ]:
Energy (eV) -7000 -6000 -5000 -4000 -3000 -2000 -1000 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 11000 12000 13000 -6000 -5500 -5000 -4500 -4000 -3500 -3000 -2500 -2000 -1500 -1000 -500 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 5500 6000 6500 7000 7500 8000 8500 9000 9500 10000 10500 11000 11500 12000 -10000 0 10000 20000 -6000 -5800 -5600 -5400 -5200 -5000 -4800 -4600 -4400 -4200 -4000 -3800 -3600 -3400 -3200 -3000 -2800 -2600 -2400 -2200 -2000 -1800 -1600 -1400 -1200 -1000 -800 -600 -400 -200 0 200 400 600 800 1000 1200 1400 1600 1800 2000 2200 2400 2600 2800 3000 3200 3400 3600 3800 4000 4200 4400 4600 4800 5000 5200 5400 5600 5800 6000 6200 6400 6600 6800 7000 7200 7400 7600 7800 8000 8200 8400 8600 8800 9000 9200 9400 9600 9800 10000 10200 10400 10600 10800 11000 11200 11400 11600 11800 12000 Si Si Ti Ti Ti Ti Ge Ge Al Al Ca Ca Ca Zn Zn Zn O O Ti Ti Ti Ca Ca Ca Ca Ge Ge Ge Zn Zn h,j,k,l,arrows,drag to pan i,o,+,-,scroll,shift-drag to zoom r,dbl-click to reset c for coordinates ? for help ? ADM-6005a_1 ADM-6005a_2 ADM-6005a_3 ADM-6005a_4 ADM-6005a_5 ADM-6005a_6 ADM-6005a_7 ADM-6005a_8 ADM-6005a_9 ADM-6005a_10 ADM-6005a_11 ADM-6005a_12 ADM-6005a_13 ADM-6005a_14 ADM-6005a_15 Spectra -100000 -80000 -60000 -40000 -20000 0 20000 40000 60000 80000 100000 120000 140000 160000 180000 -65000 -60000 -55000 -50000 -45000 -40000 -35000 -30000 -25000 -20000 -15000 -10000 -5000 0 5000 10000 15000 20000 25000 30000 35000 40000 45000 50000 55000 60000 65000 70000 75000 80000 85000 90000 95000 100000 105000 110000 115000 120000 125000 -100000 0 100000 200000 -64000 -62000 -60000 -58000 -56000 -54000 -52000 -50000 -48000 -46000 -44000 -42000 -40000 -38000 -36000 -34000 -32000 -30000 -28000 -26000 -24000 -22000 -20000 -18000 -16000 -14000 -12000 -10000 -8000 -6000 -4000 -2000 0 2000 4000 6000 8000 10000 12000 14000 16000 18000 20000 22000 24000 26000 28000 30000 32000 34000 36000 38000 40000 42000 44000 46000 48000 50000 52000 54000 56000 58000 60000 62000 64000 66000 68000 70000 72000 74000 76000 78000 80000 82000 84000 86000 88000 90000 92000 94000 96000 98000 100000 102000 104000 106000 108000 110000 112000 114000 116000 118000 120000 122000 124000 126000 Counts

To quantify this data, we need standard spectra. The standard spectra are preprocessed to optimize the fitting process. The spectrum files are read from disk and must contain the following properties :BeamEnergy, :LiveTime, :ProbeCurrent and :TakeOffAngle to be fit and quantified. It is possible to pre-load the standard spectra and edit the properties in the script if they are not present in the spectrum file.

In [ ]:
path=joinpath(datadir(), "exp_raw", "ADM6005a spectra")
refs=references( [
        reference(n"O", joinpath(path, "SiO2 std.msa"), mat"SiO2"),
        reference(n"Si", joinpath(path, "SiO2 std.msa"), mat"SiO2"),
        reference(n"Al", joinpath(path, "Al std.msa"), mat"Al"),
        reference(n"Ca", joinpath(path, "CaF2 std.msa"), mat"CaF2"),
        reference(n"Ti", joinpath(path, "Ti trimmed.msa"), mat"Ti"),
        reference(n"Zn", joinpath(path, "Zn std.msa"), mat"Zn"),
        reference(n"Ge", joinpath(path, "Ge std.msa"), mat"Ge"),
    ], 132.0
)
using DataFrames
asa(DataFrame, refs)  # Describe the references in a DataFrame
Out[ ]:
12×9 DataFrame
RowSpectrumBeam Energy (keV)Probe Current (nA)Live Time (s)MaterialLinesROIP-to-BS-to-N
StringFloat64Float64Float64Material…Array…UnitRang…Float64Float64
1SiO220.05.0100.0SiO2[Si=0.4674,O=0.5326]O K-L3 + 1 other175:22821.87884119.94
2SiO220.05.0100.0SiO2[Si=0.4674,O=0.5326]Si K-L3 + 3 others409:49120.09388705.35
3Al20.05.0100.001Al[Al=1.0000]Al K-L3 + 3 others360:43319.853714014.8
4CaF220.02.5100.0CaF2[Ca=0.5133,F=0.4867]Ca K-L3 + 3 others788:93819.44665767.9
5Ti20.05.0100.001Ti[Ti=1.0000]Ti L3-M5 + 13 others153:2264.70571092.39
6Ti20.05.0100.001Ti[Ti=1.0000]Ti K-L3 + 3 others948:112522.53310719.4
7Zn20.05.0100.002Zn[Zn=1.0000]Zn L3-M5 + 13 others247:34817.61779149.2
8Zn20.05.0100.002Zn[Zn=1.0000]Zn K-L3 + 1 other1753:188213.614272.31
9Zn20.05.0100.002Zn[Zn=1.0000]Zn K-M3 + 3 others1945:20652.6177686.219
10Ge20.05.0100.001Ge[Ge=1.0000]Ge L3-M5 + 15 others277:39714.90059088.97
11Ge20.05.0100.001Ge[Ge=1.0000]Ge K-L3 + 1 other1997:213510.57293025.67
12Ge20.05.0100.001Ge[Ge=1.0000]Ge K-M3 + 5 others2223:23581.93878478.809

The references are then used to fit the unknown spectra and matrix correction is performed. The quantify(...) method can be used to perform both operations or the fit_spectrum(...) method can be used to fit the references to the unknowns to produce k-ratios and then the k-ratios fed to the quantify(...) method to perform matrix correction. The later two-step process is useful when you want to access the k-ratio data.

In [ ]:
qrs=map(sp->quantify(sp, refs),specs)
qdf=asa(DataFrame, qrs, nominal=adm)
Out[ ]:
16×9 DataFrame
RowMaterialOAlSiCaTiZnGeTotal
StringAbstract…Abstract…Abstract…Abstract…Abstract…Abstract…Abstract…Abstract…
1ADM-6005a_10.3807590.06802230.03882770.06394990.07335950.1142640.3188521.05803
2ADM-6005a_20.3796430.06789270.03890750.06393140.07289920.1140220.3195171.05681
3ADM-6005a_30.3817520.06809040.03870230.06404760.07305050.1141830.3175111.05734
4ADM-6005a_40.380670.0680380.03895730.06379970.0730060.1141540.318971.05759
5ADM-6005a_50.3808510.06835360.03866460.06400860.07307130.1141850.3177271.05686
6ADM-6005a_60.3813980.06828510.03881310.06369240.07297450.1137070.3171951.05607
7ADM-6005a_70.3821090.06785880.03891130.06402670.07322650.1136660.3190451.05884
8ADM-6005a_80.3811850.0682170.03903830.0640360.07311370.1138160.3194571.05886
9ADM-6005a_90.3812490.06817750.03896370.06397980.07304210.1136190.3165891.05562
10ADM-6005a_100.3803730.06830980.03893390.0637950.07283260.1142850.3180741.0566
11ADM-6005a_110.3810090.06807630.03904290.06398640.07295920.1139970.3188051.05788
12ADM-6005a_120.3817740.0679320.0389310.06392850.07275760.1139030.3178321.05706
13ADM-6005a_130.3816880.06814770.03882430.06408710.07294260.1138030.3187031.0582
14ADM-6005a_140.381680.06782880.03898960.06399930.07287090.1137780.3182761.05742
15ADM-6005a_150.3823180.06829550.03879060.06404590.07310220.1142170.3186651.05943
16ADM-6005a0.33990.06640.04050.06830.07130.10550.30370.9956

We can use the describe(...) method from DataFrames to generate summary statistics.

In [ ]:
describe(qdf[1:end-1,2:end], :mean, :std, :min, :max)
Out[ ]:
8×5 DataFrame
Rowvariablemeanstdminmax
SymbolFloat64Float64Float64Float64
1O0.3812310.0007045080.3796430.382318
2Al0.06810170.0001724220.06782880.0683536
3Si0.03888660.0001132230.03866460.0390429
4Ca0.06395430.0001109570.06369240.0640871
5Ti0.07301390.0001534520.07275760.0733595
6Zn0.1139730.00023190.1136190.114285
7Ge0.3183480.0008457660.3165890.319517
8Total1.057510.001063291.055621.05943

Or we can take a closer look at the data from one spectrum.

In [ ]:
asa(DataFrame, qrs[1])
Out[ ]:
7×14 DataFrame
RowLabelElementStandardXraysCΔCkΔkgZAFcgZAFc
StringStringStringStringFloat64Float64Float64?Float64?Float64?Float64?Float64?Float64?Float64?Float64?
1ADM-6005a_1TiTiTi K-L3 + 3 others0.073360.000120.06432760.0001032821.00.9256340.9473281.01.00.876879
2ADM-6005a_1AlAlAl K-L3 + 3 others0.0680220.000120.02807664.99883e-51.01.033960.3986851.00131.00.412758
3ADM-6005a_1GeGeGe K-L3 + 1 other0.318850.000640.2634790.0005254191.00.8314310.9938681.01.00.826333
4ADM-6005a_1CaCaF2Ca K-L3 + 3 others0.063958.9e-50.1213890.0001686611.01.043120.9265011.008221.00.974395
5ADM-6005a_1ZnZnZn K-L3 + 1 other0.114260.000280.1117120.0002724371.00.8793670.9995011.112331.00.977661
6ADM-6005a_1SiSiO2Si K-L3 + 3 others0.0388289.1e-50.05305520.0001241761.01.110120.5750321.000571.00.63872
7ADM-6005a_1OSiO2O K-L3 + 1 other0.380760.000530.4875580.0006728071.01.109080.6149520.9998681.00.681942

Associated with each spectrum are a collection of properties. These properties are often read from the spectrum file but can also be assigned manually.

In [ ]:
NeXLCore.properties(sp)
Out[ ]:
Dict{Symbol, Any} with 11 entries:
  :Owner           => "Bruker Nano"
  :RealTime        => 112.108
  :BeamEnergy      => 20000.0
  :Title           => "Bruker Nano spectrum Acquisition"
  :Elevation       => 0.698132
  :LiveTime        => 100.001
  :Filename        => "c:\\Users\\nritchie\\Desktop\\IUMAS\\data\\exp_raw\\ADM6…
  :ProbeCurrent    => 5.0
  :TakeOffAngle    => 0.698132
  :Name            => "ADM-6005a_1"
  :AcquisitionTime => DateTime("2019-07-12T09:33:00")

Properties are identified with a Symbol which is represented by a colon ':' plus a name. There are many predefined spectrum properties but you can also define your own to associate custom data with spectra.

In [ ]:
sp[:BeamEnergy], sp[:ProbeCurrent], sp[:LiveTime], sp[:TakeOffAngle], dose(sp)
Out[ ]:
(20000.0, 5.0, 100.001, 0.6981317007977318, 500.005)

This is one way to index spectra to extract data from the spectrum. It is also possible to extract channel count data using various different types as indices. To demonstrate this, we will use Integer, AbstractFloat or CharXRay types to index the same channel in the spectrum - the channel associated with the Ca K-L3 peak.

In [ ]:
channel(n"Ca K-L3", sp), channel(Float64, n"Ca K-L3", sp), channel(energy(n"Ca K-L3"), sp), channel(Float64, energy(n"Ca K-L3"), sp)
Out[ ]:
(833, 833.384891053903, 833, 833.384891053903)

In the first case, we index the channel directly by channel index. In the second case, we use a floating point number which is assumed to represent an energy to index the same channel. Finally, we use a CharXRay to index the same channel.

In [ ]:
sp[833], sp[energy(n"Ca K-L3")], sp[n"Ca K-L3"]
Out[ ]:
(18042.0, 18042.0, 18042.0)

Often it is useful to be able to determine which of a collection of spectra collected from a single material under similar conditions are similar to one another. The notion being that the 'best' spectra are those that are most similar. These can then be summed together to produce a single spectrum that best represents the material. The findsimilar(...) method extracts the similar spectra from a list and removes outliers.

In [ ]:
plot(sum(findsimilar(specs, refs.detector, n"O")), klms=[n"O", n"Al", n"Si", n"Ca", n"Ti", n"Zn", n"Ge"])
Out[ ]:
Energy (eV) -25000 -20000 -15000 -10000 -5000 0 5000 10000 15000 20000 25000 30000 35000 40000 45000 -20000 -18000 -16000 -14000 -12000 -10000 -8000 -6000 -4000 -2000 0 2000 4000 6000 8000 10000 12000 14000 16000 18000 20000 22000 24000 26000 28000 30000 32000 34000 36000 38000 40000 -20000 0 20000 40000 -20000 -19000 -18000 -17000 -16000 -15000 -14000 -13000 -12000 -11000 -10000 -9000 -8000 -7000 -6000 -5000 -4000 -3000 -2000 -1000 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 11000 12000 13000 14000 15000 16000 17000 18000 19000 20000 21000 22000 23000 24000 25000 26000 27000 28000 29000 30000 31000 32000 33000 34000 35000 36000 37000 38000 39000 40000 Si Si Ti Ti Ti Ti Ge Ge Al Al Ca Ca Ca Zn Zn Zn O O Ti Ti Ti Ca Ca Ca Ca Ge Ge Ge Zn Zn h,j,k,l,arrows,drag to pan i,o,+,-,scroll,shift-drag to zoom r,dbl-click to reset c for coordinates ? for help ? Sum[15 spectra] Spectrum -1500000 -1000000 -500000 0 500000 1000000 1500000 2000000 2500000 -1000000 -900000 -800000 -700000 -600000 -500000 -400000 -300000 -200000 -100000 0 100000 200000 300000 400000 500000 600000 700000 800000 900000 1000000 1100000 1200000 1300000 1400000 1500000 1600000 1700000 1800000 1900000 -1000000 0 1000000 2000000 -950000 -900000 -850000 -800000 -750000 -700000 -650000 -600000 -550000 -500000 -450000 -400000 -350000 -300000 -250000 -200000 -150000 -100000 -50000 0 50000 100000 150000 200000 250000 300000 350000 400000 450000 500000 550000 600000 650000 700000 750000 800000 850000 900000 950000 1000000 1050000 1100000 1150000 1200000 1250000 1300000 1350000 1400000 1450000 1500000 1550000 1600000 1650000 1700000 1750000 1800000 1850000 Counts

The similarity is determined by comparing each spectrum with the sum of the other spectra and calculating a score that respects count statistics. Similar spectra have a similarity(...) score of approximately unity.

In [ ]:
DataFrame(Name=name.(specs), S1=NeXLSpectrum.similarity(specs, refs.detector, adm), S2=NeXLSpectrum.similarity(specs), Counts=integrate.(specs))
Out[ ]:
15×4 DataFrame
RowNameS1S2Counts
StringFloat64Float64Float64
1ADM-6005a_10.8614811.041116.81189e6
2ADM-6005a_20.8520620.9704976.81126e6
3ADM-6005a_30.8546750.999766.81006e6
4ADM-6005a_40.8542780.9802646.81625e6
5ADM-6005a_50.890410.9806026.80921e6
6ADM-6005a_60.8909280.9624676.8086e6
7ADM-6005a_70.8412130.9766286.81929e6
8ADM-6005a_80.7788260.9711666.81616e6
9ADM-6005a_90.838671.016526.816e6
10ADM-6005a_100.8473911.003536.81184e6
11ADM-6005a_110.8616670.9754456.81146e6
12ADM-6005a_120.8231051.014216.81786e6
13ADM-6005a_130.8143980.9742836.81113e6
14ADM-6005a_140.7775660.9723486.81502e6
15ADM-6005a_150.7807070.9424426.81972e6

HyperSpectra with NeXLSpectrum¶

A HyperSpectrum represents a multi-dimensional array of point spectra. NeXLSpectrum can handle a arbitrary number of dimensions like 1-D (linescans), 2-D (area scans) or higher dimensions.

There are currently three ways to load a hyperspectrum.

  • From a Lispix-style RPL/RAW file pair
  • From a SEMantics-style PTZ file
  • From a HyperSpy-style HDF5 file

We will use the DataDeps library to download the hyperspectrum data on demand from the NIST data website. You will need to be connected to the Internet.

In [ ]:
# DataDeps is a utility for downloading data on demand from a URL
using DataDeps

# Where do I want to put the data?
ENV["DATADEPS_LOAD_PATH"] = joinpath(datadir(), "exp_raw")
ENV["DATADEPS_ALWAYS_ACCEPT"]="true"

# Register the data sets using the names "MnDeepSeaNodule" and "MnDeepSeaNoduleStandards"
register(DataDep("MnDeepSeaNodule",
    """
    Dataset:       Deep Sea Manganese Nodule
    Acquired by:   Nicholas W.M. Ritchie
    License:       CC-SA 3.0
    """,
    raw"https://data.nist.gov/od/ds/mds2-2467/MnNodule.tar.gz",
    "5b5b6623b8f4daca3ff3073708442ac5702ff690aa12668659875ec5642b458d",
    post_fetch_method=DataDeps.unpack
));

register(DataDep("MnDeepSeaNoduleStandards",
    """
    Dataset:       Standard Spectra for Deep Sea Manganese Nodule Dataset
    Acquired by:   Nicholas W.M. Ritchie
    License:       CC-SA 3.0
    """,
    raw"https://data.nist.gov/od/ds/mds2-2467/MnNodule_Standards.tar.gz",
    "69283ba72146932ba451e679cf02fbd6b350f96f6d012d50f589ed9dd2e35f1a",
    post_fetch_method=DataDeps.unpack
));

The data is in RPL/RAW format - a format designed by David Bright at NIST and available from many EDS vendor's software. In addition to the data in the RPL/RAW format, we also need to provide spectrum meta-data and the detector energy scale.

In [ ]:
raw = readrplraw(joinpath(datadep"MnDeepSeaNodule","map[15]"))
les = LinearEnergyScale(0.0, 10.0)
props = Dict{Symbol,Any}(
    :LiveTime => 0.72*4.0*18.0*3600.0/(1024.0*1024.0),
    :BeamEnergy => 20.0e3,
    :TakeOffAngle => deg2rad(35.0),
    :ProbeCurrent => 1.0,
    :Name => "Deep Sea Mn Nodule",
)
hs = HyperSpectrum(les, props, raw, fov = (0.2, 0.2))
Out[ ]:

1024 × 1024 HyperSpectrum{UInt16,(:Y, :X)}[Deep Sea Mn Nodule), 0.0 + 10.0⋅ch eV, 2048 ch]

It is easy to generate crude elemental maps from HyperSpectrum structs. These are simple ROI maps constructed by summing a range of channels around the central channel for the specified X-ray transition.

In [ ]:
hs[n"Mn K-L3"]
Out[ ]:

Or I can create a RGB-colorized image where R represents an ROI around the first chararacteristic X-ray, green the second and blue the third.

In [ ]:
hs[n"Mn K-L3", n"Si K-L3", n"O K-L3"]
Out[ ]:

Much like spectra, hyperspectra have properties. However, the same set of properties apply to all spectra within a hyperspectrum. There is one exception, :LiveTime, which can be an array of the same dimension as the hyperspectrum to represent a unique live-time per pixel.

In [ ]:
NeXLCore.properties(hs)
Out[ ]:
Dict{Symbol, Any} with 5 entries:
  :Name         => "Deep Sea Mn Nodule"
  :BeamEnergy   => 20000.0
  :LiveTime     => 0.177979
  :TakeOffAngle => 0.610865
  :ProbeCurrent => 1.0

To speed up processing bin the data using a block size > 1. 4 will provide 16x faster processing but lower resolution results.

In [ ]:
# This speeds up the processing by binning a 1024 x 1024 hyperspectrum to 1024/block_size x 1024/block_size.
block_size = 1
hs= block_size > 1 ? block(hs, (block_size, block_size)) : hs
Out[ ]:

1024 × 1024 HyperSpectrum{UInt16,(:Y, :X)}[Deep Sea Mn Nodule), 0.0 + 10.0⋅ch eV, 2048 ch]

Unfortunately, ROI maps suffer from many shortcomings including

  • not being background corrected to eliminate the continuum signal
  • not being scaled for differences in generation efficiency between different characteristic X-rays.

While they are convenient as a first glimse at the data, they can be very deceptive and should rarely be included in reports.

Instead, it is better to fit the point spectra within the hyperspectrum using measured references.

In [ ]:
hs_elms = [ n"Ag", n"Al", n"Ba", n"C", n"Ca", n"Ce", n"Cr", n"Cu", n"Fe", n"S", n"P", n"K", n"Mg", n"O", n"Mn", n"Na", n"Cl", n"Ni", n"Si", n"Ti", n"Zn" ]
plot(maxpixel(hs), klms = hs_elms, xmax=10.0e3)
Out[ ]:
Energy (eV) -15000 -10000 -5000 0 5000 10000 15000 20000 25000 -10000 -9000 -8000 -7000 -6000 -5000 -4000 -3000 -2000 -1000 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 11000 12000 13000 14000 15000 16000 17000 18000 19000 20000 -10000 0 10000 20000 -10000 -9500 -9000 -8500 -8000 -7500 -7000 -6500 -6000 -5500 -5000 -4500 -4000 -3500 -3000 -2500 -2000 -1500 -1000 -500 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 5500 6000 6500 7000 7500 8000 8500 9000 9500 10000 10500 11000 11500 12000 12500 13000 13500 14000 14500 15000 15500 16000 16500 17000 17500 18000 18500 19000 19500 20000 Si Si Ti Ti Ti K K K K Ti Ti Ti Ti Ag Ag Ag Ag Zn Zn Cu Cu Cu Ag Ag Ag Ag Ag Mn Mn Mn Cu Cu Al Al P P Mn Mn Mn Fe Fe Fe S S Ca Ca Ca Fe Fe Fe Mg Mg Cr Cr Cr Ni Ni Ni O O Cl Cl Zn Zn Zn Ba Ba Ba Ba Ba Ca Ca Ca Ca C C Ce Ce Ce Ce Ce Ba Ba Ba Ba K K Na Na Cr Cr Cr Ce Ce Ce Ce Ni Ni h,j,k,l,arrows,drag to pan i,o,+,-,scroll,shift-drag to zoom r,dbl-click to reset c for coordinates ? for help ? MaxPixel[Deep Sea Mn Nodule] Spectrum -700 -600 -500 -400 -300 -200 -100 0 100 200 300 400 500 600 700 800 900 1000 1100 1200 1300 -550 -500 -450 -400 -350 -300 -250 -200 -150 -100 -50 0 50 100 150 200 250 300 350 400 450 500 550 600 650 700 750 800 850 900 950 1000 1050 -1000 0 1000 2000 -520 -500 -480 -460 -440 -420 -400 -380 -360 -340 -320 -300 -280 -260 -240 -220 -200 -180 -160 -140 -120 -100 -80 -60 -40 -20 0 20 40 60 80 100 120 140 160 180 200 220 240 260 280 300 320 340 360 380 400 420 440 460 480 500 520 540 560 580 600 620 640 660 680 700 720 740 760 780 800 820 840 860 880 900 920 940 960 980 1000 1020 Counts
In [ ]:
print(symbol.(sort(hs_elms)))
["C", "O", "Na", "Mg", "Al", "Si", "P", "S", "Cl", "K", "Ca", "Ti", "Cr", "Mn", "Fe", "Ni", "Cu", "Zn", "Ag", "Ba", "Ce"]
In [ ]:
refs = references( [
        reference(n"C", joinpath(datadep"MnDeepSeaNoduleStandards", "C std.msa")),
        reference(n"O", joinpath(datadep"MnDeepSeaNoduleStandards", "MgO std.msa")),
        reference(n"Na", joinpath(datadep"MnDeepSeaNoduleStandards", "NaCl std.msa")),
        reference(n"Mg", joinpath(datadep"MnDeepSeaNoduleStandards", "Mg std.msa")),
        reference(n"Al", joinpath(datadep"MnDeepSeaNoduleStandards", "Al std.msa")),
        reference(n"Si", joinpath(datadep"MnDeepSeaNoduleStandards", "Si std.msa")),
        reference(n"P", joinpath(datadep"MnDeepSeaNoduleStandards", "GaP std.msa")), 
        reference(n"S", joinpath(datadep"MnDeepSeaNoduleStandards", "FeS2 std.msa")), 
        reference(n"Cl", joinpath(datadep"MnDeepSeaNoduleStandards", "NaCl std.msa")), 
        reference(n"K", joinpath(datadep"MnDeepSeaNoduleStandards", "KBr std.msa")), 
        reference(n"Ca", joinpath(datadep"MnDeepSeaNoduleStandards", "CaF2 std.msa")), 
        reference(n"Ti", joinpath(datadep"MnDeepSeaNoduleStandards", "Ti std.msa")), 
        reference(n"Cr", joinpath(datadep"MnDeepSeaNoduleStandards", "Cr std.msa")), 
        reference(n"Mn", joinpath(datadep"MnDeepSeaNoduleStandards", "Mn std.msa")), 
        reference(n"Fe", joinpath(datadep"MnDeepSeaNoduleStandards", "Fe std.msa")), 
        reference(n"Ni", joinpath(datadep"MnDeepSeaNoduleStandards", "Ni std.msa")), 
        reference(n"Cu", joinpath(datadep"MnDeepSeaNoduleStandards", "Cu std.msa")), 
        reference(n"Zn", joinpath(datadep"MnDeepSeaNoduleStandards", "Zn std.msa")), 
        reference(n"Ag", joinpath(datadep"MnDeepSeaNoduleStandards", "Ag std.msa")), 
        reference(n"Ba", joinpath(datadep"MnDeepSeaNoduleStandards", "BaF2 std.msa"))     
    ], 132.0)

asa(DataFrame, refs)
[ Info: A material containing Ba and F is not a suitable reference for "Ba M5-N3 + 16 others" due to a peak interference.
Out[ ]:
33×9 DataFrame
8 rows omitted
RowSpectrumBeam Energy (keV)Probe Current (nA)Live Time (s)MaterialLinesROIP-to-BS-to-N
SubStrin…Float64Float64Float64Material…Array…UnitRang…Float64Float64
1C std20.00.973721184.6C[C=1.0000,1.90 g/cm³]C K-L3 + 1 other16:4024.14315279.7
2MgO std20.00.97481161.48MgO[Mg=0.6030,O=0.3970,5.00 g/cm³]O K-L3 + 1 other39:6612.78898706.54
3NaCl std20.00.974541156.51NaCl[Na=0.3934,Cl=0.6066,5.00 g/cm³]Na K-L3 + 1 other89:1206.218399555.21
4Mg std20.00.974751149.05Mg[Mg=1.0000,1.74 g/cm³]Mg K-L3 + 1 other110:14217.447134147.0
5Al std20.00.973961151.25Al[Al=1.0000,1.00 g/cm³]Al K-L3 + 3 others132:16919.638936840.4
6Si std20.00.974651149.58Si[Si=1.0000,2.95 g/cm³]Si K-L3 + 3 others157:19825.039340461.5
7GaP std20.00.974511161.35GaP[P=0.3076,Ga=0.6924,5.00 g/cm³]P K-L3 + 3 others184:2308.271729393.98
8FeS2 std20.00.973341160.73FeS2[S=0.5345,Fe=0.4655,5.00 g/cm³]S K-L3 + 3 others212:26416.838821532.6
9NaCl std20.00.974541156.51NaCl[Na=0.3934,Cl=0.6066,5.00 g/cm³]Cl K-L3 + 3 others243:30023.51225791.1
10KBr std20.00.973921147.26KBr[Br=0.6714,K=0.3286,5.00 g/cm³]K K-L3 + 3 others310:3797.045759501.98
11CaF2 std20.00.974011167.54CaF2[Ca=0.5133,F=0.4867,3.18 g/cm³]Ca K-L3 + 3 others347:42217.527419125.9
12Ti std20.00.97521160.26Ti[Ti=1.0000,5.00 g/cm³]Ti L3-M5 + 13 others28:654.350943026.39
13Ti std20.00.97521160.26Ti[Ti=1.0000,5.00 g/cm³]Ti K-L3 + 3 others427:51620.563626561.6
⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮
22Fe std20.00.974361164.66Fe[Fe=1.0000,7.87 g/cm³]Fe K-M3 + 3 others680:7323.778543010.53
23Ni std20.00.972581164.08Ni[Ni=1.0000,5.00 g/cm³]Ni L3-M5 + 13 others62:10810.005211386.8
24Ni std20.00.972581164.08Ni[Ni=1.0000,5.00 g/cm³]Ni K-L3 + 1 other718:77816.041314413.1
25Ni std20.00.972581164.08Ni[Ni=1.0000,5.00 g/cm³]Ni K-M3 + 3 others799:8553.183282401.56
26Cu std20.00.972421159.85Cu[Cu=1.0000,8.90 g/cm³]Cu L3-M5 + 13 others69:11715.370319251.1
27Cu std20.00.972421159.85Cu[Cu=1.0000,8.90 g/cm³]Cu K-L3 + 1 other773:83614.093612296.1
28Cu std20.00.972421159.85Cu[Cu=1.0000,8.90 g/cm³]Cu K-M3 + 3 others862:9202.930672123.89
29Zn std20.00.972921157.03Zn[Zn=1.0000,5.00 g/cm³]Zn L3-M5 + 13 others76:12617.349723748.8
30Zn std20.00.972921157.03Zn[Zn=1.0000,5.00 g/cm³]Zn K-L3 + 1 other831:89612.506910519.0
31Zn std20.00.972921157.03Zn[Zn=1.0000,5.00 g/cm³]Zn K-M3 + 3 others928:9882.675581862.94
32Ag std20.00.98647921.028Ag[Ag=1.0000,5.00 g/cm³]Ag L3-M5 + 25 others247:3935.2894511615.2
33BaF2 std20.00.989861171.3BaF2[Ba=0.7833,F=0.2167,5.00 g/cm³]Ba L3-M5 + 29 others377:6172.854567169.28

Now we will fit the references to the hyperspectrum unknown. This is simple enough:

krs = fit_spectrum(hs, refs, mode = :Fast)

This is time consuming so once we've done it, we want to store the result in a file and reload the result rather than recalculate it. DrWatsons provides a mechanism to do this based on a function produce_or_load(....). You specify where to put the data and it builds a filename based on the prefix, suffix and the configuration paramters in @dict(mode, x_dim, y_dim). The suffix specifies that the data is written using the JLD2 library in HDF5 format.

In [ ]:
mode = :Full  # Single threaded took 92 seconds for :Fast or 6000 seconds for :Full (many cores take about 1/6 this.)

# krs = @time fit_spectrum(hs, refs, mode = mode)  

using FileIO, Parameters
x_dim, y_dim = size(hs)
data, file = produce_or_load(
        joinpath(datadir(),"exp_pro"), # Path
        @dict(mode, x_dim, y_dim),     # Configuration
        prefix="kratios", suffix="jld2" # Filename parameters
    ) do config
        Dict("kratios" => @time fit_spectrum(hs, refs, mode = mode))
end
@show file
krs = data["kratios"]
file = "c:\\Users\\nritchie\\Desktop\\IUMAS\\data\\exp_pro\\kratios_mode=Full_x_dim=1024_y_dim=1024.jld2"
Out[ ]:
33-element Vector{KRatios}:
 k[C, C K-L3 + 1 other] = Float32[ (1024, 1024) ]
 k[MgO, O K-L3 + 1 other] = Float32[ (1024, 1024) ]
 k[NaCl, Na K-L3 + 1 other] = Float32[ (1024, 1024) ]
 k[Mg, Mg K-L3 + 1 other] = Float32[ (1024, 1024) ]
 k[Al, Al K-L3 + 3 others] = Float32[ (1024, 1024) ]
 k[Si, Si K-L3 + 3 others] = Float32[ (1024, 1024) ]
 k[GaP, P K-L3 + 3 others] = Float32[ (1024, 1024) ]
 k[FeS2, S K-L3 + 3 others] = Float32[ (1024, 1024) ]
 k[NaCl, Cl K-L3 + 3 others] = Float32[ (1024, 1024) ]
 k[KBr, K K-L3 + 3 others] = Float32[ (1024, 1024) ]
 k[CaF2, Ca K-L3 + 3 others] = Float32[ (1024, 1024) ]
 k[Ti, Ti L3-M5 + 13 others] = Float32[ (1024, 1024) ]
 k[Ti, Ti K-L3 + 3 others] = Float32[ (1024, 1024) ]
 ⋮
 k[Fe, Fe K-M3 + 3 others] = Float32[ (1024, 1024) ]
 k[Ni, Ni L3-M5 + 13 others] = Float32[ (1024, 1024) ]
 k[Ni, Ni K-L3 + 1 other] = Float32[ (1024, 1024) ]
 k[Ni, Ni K-M3 + 3 others] = Float32[ (1024, 1024) ]
 k[Cu, Cu L3-M5 + 13 others] = Float32[ (1024, 1024) ]
 k[Cu, Cu K-L3 + 1 other] = Float32[ (1024, 1024) ]
 k[Cu, Cu K-M3 + 3 others] = Float32[ (1024, 1024) ]
 k[Zn, Zn L3-M5 + 13 others] = Float32[ (1024, 1024) ]
 k[Zn, Zn K-L3 + 1 other] = Float32[ (1024, 1024) ]
 k[Zn, Zn K-M3 + 3 others] = Float32[ (1024, 1024) ]
 k[Ag, Ag L3-M5 + 25 others] = Float32[ (1024, 1024) ]
 k[BaF2, Ba L3-M5 + 29 others] = Float32[ (1024, 1024) ]

The result of fit_spectrum(...) is a set k-ratio maps. Which can be readily visualized and represent a more quantitative perspective on the spectrum data. The principle advantage being that the data is background corrected so that continuum is not mistaken for a small quantity of an element.

To visualize the k-ratios, it is best to pick one k-ratio per element using optimizeks(...) and then normalize the k-ratios on a point-by-point basis using normalizek(...). This ensures that all the k-ratios are between 0 and 1. The function asa(ElementalMap, ...) facilitates this process by converting a Vector{KRatios} into a dictionary of images indexed by Element instances.

In [ ]:
using Colors
imgs = asa(ElementalMap, krs, Gray)
imgs[n"Mn"]
Out[ ]:

We can perform matrix correction on the Vector{KRatios} using the quantify(...) method. It returns a Materials struct which is a memory efficient way of representing a multi-dimensional array of Material.

In [ ]:
# quant=quantify(krs, name=hs[:Name], ty=Float32) # Takes ~80 minute for 1024 x 1024 (~5 ms per pixel)

data, file = produce_or_load(
        joinpath(datadir(),"exp_pro"), 
        @dict(mode, x_dim, y_dim),
        prefix="quantify", suffix="jld2"
    ) do config
        Dict("mass_fractions" => @time quantify(krs, name=hs[:Name], ty=Float32))
end
@show file
quant = data["mass_fractions"]
file = "c:\\Users\\nritchie\\Desktop\\IUMAS\\data\\exp_pro\\quantify_mode=Full_x_dim=1024_y_dim=1024.jld2"
Out[ ]:
Materials[Deep Sea Mn Nodule, 1024 × 1024 × (C, O, Na, Mg, Al, Si, P, S, Cl, K, Ca, Ti, Cr, Mn, Fe, Ni, Cu, Zn, Ag, Ba)]

The Materials struct can be indexed at a specific coordinate to extract the composition at that coordinate as a Material.

In [ ]:
quant[64,192]
Out[ ]:
Deep Sea Mn Nodule[64,192][S=0.0058,Ca=0.0172,Mg=0.0224,Zn=0.0052,Ni=0.0232,O=0.4625,Cl=0.0008,K=0.0067,C=0.1091,Na=0.0150,Si=0.0231,Cu=0.0136,Mn=0.3366,Al=0.0215,Fe=0.0336]

Or we can index the results using an Element in which case, we pull out a slice representing the mass fractions associated with that element.

In [ ]:
t=quant[n"Mn"]
Out[ ]:
1024×1024 Matrix{Float32}:
 0.223394    0.311174   0.332411    …  0.192086    0.18319     0.182768
 0.269441    0.277169   0.307894       0.187049    0.187401    0.198357
 0.27579     0.272186   0.242913       0.1994      0.215164    0.195631
 0.284984    0.254272   0.251819       0.221909    0.19468     0.228718
 0.271912    0.248145   0.231546       0.181427    0.315283    0.155952
 0.00486168  0.235095   0.232646    …  0.209557    0.0161902   0.198351
 0.216475    0.192771   0.306901       0.200423    0.165684    0.218023
 0.220113    0.231561   0.285585       0.213188    0.22912     0.21971
 0.248034    0.252878   0.272362       0.19302     0.214149    0.281862
 0.24977     0.260708   0.239762       0.301234    0.185244    0.209922
 0.266883    0.271193   0.297616    …  0.165329    0.190385    0.211984
 0.301004    0.146325   0.274044       0.26669     0.218611    0.207211
 0.29647     0.303459   0.285447       0.265315    0.248976    0.217168
 ⋮                                  ⋱                          
 0.324508    0.324034   0.26046        0.00610759  0.00819175  0.0323414
 0.223967    0.298605   0.231759       0.247855    0.0136075   0.0366124
 0.267053    0.2395     0.00541237     0.00076621  0.0255185   0.0395083
 0.0377666   0.161325   0.092179    …  0.0         0.416347    0.00254732
 0.147561    0.097022   0.0636286      0.193515    0.00693751  0.00533547
 0.0941413   0.0659752  0.0677662      0.00869319  0.00985284  0.00891134
 0.101526    0.07201    0.104355       0.0264528   0.017331    0.00477545
 0.27833     0.129675   0.188148       0.036561    0.0274113   0.0151129
 0.15678     0.207637   0.274506    …  0.0463395   0.0424375   0.334446
 0.225777    0.146953   0.385326       0.069456    0.0334889   0.0160826
 0.208606    0.294549   0.346879       0.0603972   0.0310457   0.0142638
 0.220558    0.30014    0.321251       0.0684462   0.0361965   0.0130948

Using broadcast, we can apply Material functions to Materials objects. The results is an Array with the same shape but with the result contents.

In [ ]:
anal_tot=analyticaltotal.(quant)
Out[ ]:
1024×1024 Matrix{Float64}:
 0.605907  0.941084  0.991103  1.10609   …  0.822585  0.779711  0.747124
 0.856864  0.925094  1.00678   0.844349     0.807512  0.747789  0.815007
 0.928448  0.981215  0.925177  0.879346     0.772628  0.79066   0.811633
 0.922476  0.847506  0.866738  0.929373     0.672841  0.705843  0.73923
 0.828224  0.829324  0.778171  0.926691     0.591353  0.887882  0.813112
 1.14474   0.882673  0.900289  0.847674  …  0.654242  0.402655  0.661255
 0.949938  0.979103  1.12348   0.920959     0.809175  0.846791  0.691027
 0.886173  0.881917  0.966901  0.937227     0.791499  0.743163  0.660464
 0.832243  0.872016  0.92212   0.938177     0.763454  0.78746   0.69703
 0.893395  0.928652  0.976265  1.06764      0.833232  0.806506  0.716307
 0.962367  0.997508  0.95539   0.950264  …  0.856741  0.817339  0.741525
 0.942974  0.952614  0.892325  0.950321     1.09738   0.796921  0.807034
 0.88609   0.874692  0.821842  0.815836     0.864204  0.871525  0.819883
 ⋮                                       ⋱                      
 0.876421  0.912032  0.924582  0.97681      0.964411  0.84925   0.799963
 0.928951  0.985055  0.937457  0.965731     0.994776  0.881218  0.897845
 1.03631   1.05844   1.15028   1.00834      0.837187  0.955848  0.946106
 1.25698   0.961709  0.852426  1.18083   …  0.915044  1.11787   0.927067
 1.05104   1.05405   1.02976   1.0058       1.26751   1.05273   1.09016
 1.0333    1.03098   0.992141  0.958689     0.828392  0.883704  1.07129
 1.11179   1.10899   1.04152   0.972553     0.602862  0.738939  0.840986
 0.959889  1.03677   0.973237  1.03496      0.747029  0.504262  0.707248
 1.00226   1.11105   0.998813  1.84589   …  0.522819  0.560106  1.09384
 1.08641   0.958066  1.16998   1.01273      0.618802  0.740887  0.947171
 0.913257  0.957598  0.952383  0.924942     0.738241  0.926618  0.963841
 0.961834  1.04834   0.996968  0.974231     0.7779    0.896684  0.932103
In [ ]:
using Statistics
mean(anal_tot), std(anal_tot)
Out[ ]:
(1.0397453436659276, 0.1712837655499593)

To convert the mass-fractions to images, we need to ensure all mass-fractions are on the range [0,1]. We can do this by normalizing the individual points to an analytical total of unity.

In [ ]:
norm_quant=asnormalized(quant)
Out[ ]:
Materials[N[Deep Sea Mn Nodule], 1024 × 1024 × (C, O, Na, Mg, Al, Si, P, S, Cl, K, Ca, Ti, Cr, Mn, Fe, Ni, Cu, Zn, Ag, Ba)]

Now let's visualize this. In Julia, it is trivial to convert matrices representing values on the range [0,1] to images.

We will extract the plane of data corresponding to an element by indexing it with an Element. Log3Band is a function that converts numbers between 0 and 1 to RGB values. The mapping is displayed below in the legend.

In [ ]:
Log3Band.(norm_quant[n"Ni"])
Out[ ]:
In [ ]:
loadlegend("Log3BandBright.png")
Out[ ]:

The Ni elemental map can be seen to represent primarily values ranging from a few percent down to zero. Albeit, this spectrum image took 18 hours to collect but the depth of information might surprise many who are used to thinking of spectrum images as a crude tool.

Let's output various perspectives from each element to the 'plots' directory. The Gray function converts values on the range [0,1] to gray-scale values.

In [ ]:
# Fast k-ratio maps
imgs = asa(ElementalMap, krs, Log3Band) 
for (el, img) in imgs
    save(joinpath(plotsdir(),"$(symbol(el)) - $mode Log3Band.png"), img)
end
# Full k-ratio maps
for (el, img) in asa(ElementalMap, krs, Gray) 
    save(joinpath(plotsdir(),"$(symbol(el)) - $mode Gray.png"), img)
end
# Full quantitative maps - Fast fit
for el in keys(norm_quant)
    pl = norm_quant[el]
    save(joinpath(plotsdir(),"$(symbol(el)) - Quant $mode Gray.png"), Gray.(pl))
    save(joinpath(plotsdir(),"$(symbol(el)) - Quant $mode Log3Band.png"), Log3Band.(pl))
end

Let's compare the results from the full fit with the results from the fast fit.

In [ ]:
els = sort([ keys(imgs)...])
print(symbol.(els))
["C", "O", "Na", "Mg", "Al", "Si", "P", "S", "Cl", "K", "Ca", "Ti", "Cr", "Mn", "Fe", "Ni", "Cu", "Zn", "Ag", "Ba"]

Plot the k-ratio images.

In [ ]:
map( el->imgs[el], els)
Out[ ]:
(a vector displayed as a row to save space)

Plot the quantitative maps.

In [ ]:
map(el->Log3Band.(norm_quant[el]), els)
Out[ ]:
(a vector displayed as a row to save space)

Clearly, the light element maps are where we see the influence of matrix correction. When X-ray absorption is small, the k-ratio is generally a good approximation of the mass-fraction. However, for low energy X-rays, the absorption can be strong leading to a significant underestimation of the mass fraction. We see this particularly in the carbon data.

In [ ]:
using Statistics
mean(analyticaltotal.(quant)), mean(analyticaltotal.(quant))
Out[ ]:
(1.0397453436659276, 1.0397453436659276)

To explore the data, let's plot histograms of the mass fractions from each element.

In [ ]:
colors = distinguishable_colors(length(els), colorant"light blue", lchoices=range(50, stop=100, length=15))
    
plot(
    ( layer(x=quant[el], Geom.histogram, Theme(default_color=colors[i]), alpha=[0.6]) for (i, el) in enumerate(els) )...,
    Guide.xlabel("Mass Fraction"), Guide.ylabel("Counts"), 
    Guide.manual_color_key("Element", symbol.(els), colors),
    Coord.cartesian(xmin=0.0, xmax=1.0, ymax=16.0e4/(block_size^2))
)
Out[ ]:
Mass Fraction -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 -1 0 1 2 -1.00 -0.95 -0.90 -0.85 -0.80 -0.75 -0.70 -0.65 -0.60 -0.55 -0.50 -0.45 -0.40 -0.35 -0.30 -0.25 -0.20 -0.15 -0.10 -0.05 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00 1.05 1.10 1.15 1.20 1.25 1.30 1.35 1.40 1.45 1.50 1.55 1.60 1.65 1.70 1.75 1.80 1.85 1.90 1.95 2.00 Ca Ti Cr Mn Fe Ni Cu Zn Ag Ba C O Na Mg Al Si P S Cl K Element h,j,k,l,arrows,drag to pan i,o,+,-,scroll,shift-drag to zoom r,dbl-click to reset c for coordinates ? for help ? -1.80×10⁵ -1.60×10⁵ -1.40×10⁵ -1.20×10⁵ -1.00×10⁵ -8.00×10⁴ -6.00×10⁴ -4.00×10⁴ -2.00×10⁴ 0 2.00×10⁴ 4.00×10⁴ 6.00×10⁴ 8.00×10⁴ 1.00×10⁵ 1.20×10⁵ 1.40×10⁵ 1.60×10⁵ 1.80×10⁵ 2.00×10⁵ 2.20×10⁵ 2.40×10⁵ 2.60×10⁵ 2.80×10⁵ 3.00×10⁵ 3.20×10⁵ 3.40×10⁵ -1.60×10⁵ -1.50×10⁵ -1.40×10⁵ -1.30×10⁵ -1.20×10⁵ -1.10×10⁵ -1.00×10⁵ -9.00×10⁴ -8.00×10⁴ -7.00×10⁴ -6.00×10⁴ -5.00×10⁴ -4.00×10⁴ -3.00×10⁴ -2.00×10⁴ -1.00×10⁴ 0 1.00×10⁴ 2.00×10⁴ 3.00×10⁴ 4.00×10⁴ 5.00×10⁴ 6.00×10⁴ 7.00×10⁴ 8.00×10⁴ 9.00×10⁴ 1.00×10⁵ 1.10×10⁵ 1.20×10⁵ 1.30×10⁵ 1.40×10⁵ 1.50×10⁵ 1.60×10⁵ 1.70×10⁵ 1.80×10⁵ 1.90×10⁵ 2.00×10⁵ 2.10×10⁵ 2.20×10⁵ 2.30×10⁵ 2.40×10⁵ 2.50×10⁵ 2.60×10⁵ 2.70×10⁵ 2.80×10⁵ 2.90×10⁵ 3.00×10⁵ 3.10×10⁵ 3.20×10⁵ -2.0×10⁵ 0 2.0×10⁵ 4.0×10⁵ -1.600×10⁵ -1.550×10⁵ -1.500×10⁵ -1.450×10⁵ -1.400×10⁵ -1.350×10⁵ -1.300×10⁵ -1.250×10⁵ -1.200×10⁵ -1.150×10⁵ -1.100×10⁵ -1.050×10⁵ -1.000×10⁵ -9.500×10⁴ -9.000×10⁴ -8.500×10⁴ -8.000×10⁴ -7.500×10⁴ -7.000×10⁴ -6.500×10⁴ -6.000×10⁴ -5.500×10⁴ -5.000×10⁴ -4.500×10⁴ -4.000×10⁴ -3.500×10⁴ -3.000×10⁴ -2.500×10⁴ -2.000×10⁴ -1.500×10⁴ -1.000×10⁴ -5.000×10³ 0 5.000×10³ 1.000×10⁴ 1.500×10⁴ 2.000×10⁴ 2.500×10⁴ 3.000×10⁴ 3.500×10⁴ 4.000×10⁴ 4.500×10⁴ 5.000×10⁴ 5.500×10⁴ 6.000×10⁴ 6.500×10⁴ 7.000×10⁴ 7.500×10⁴ 8.000×10⁴ 8.500×10⁴ 9.000×10⁴ 9.500×10⁴ 1.000×10⁵ 1.050×10⁵ 1.100×10⁵ 1.150×10⁵ 1.200×10⁵ 1.250×10⁵ 1.300×10⁵ 1.350×10⁵ 1.400×10⁵ 1.450×10⁵ 1.500×10⁵ 1.550×10⁵ 1.600×10⁵ 1.650×10⁵ 1.700×10⁵ 1.750×10⁵ 1.800×10⁵ 1.850×10⁵ 1.900×10⁵ 1.950×10⁵ 2.000×10⁵ 2.050×10⁵ 2.100×10⁵ 2.150×10⁵ 2.200×10⁵ 2.250×10⁵ 2.300×10⁵ 2.350×10⁵ 2.400×10⁵ 2.450×10⁵ 2.500×10⁵ 2.550×10⁵ 2.600×10⁵ 2.650×10⁵ 2.700×10⁵ 2.750×10⁵ 2.800×10⁵ 2.850×10⁵ 2.900×10⁵ 2.950×10⁵ 3.000×10⁵ 3.050×10⁵ 3.100×10⁵ 3.150×10⁵ 3.200×10⁵ Counts
In [ ]:
plot(
    ( layer(x=quant[el], Geom.histogram, Theme(default_color=colors[i]), alpha=[0.6]) for (i, el) in enumerate(els) )...,
    Guide.xlabel("Mass Fraction"), Guide.ylabel("Counts"), 
    Guide.manual_color_key("Element", symbol.(els), colors),
    Coord.cartesian(xmin=0.0, xmax=0.05, ymax=16.0e4/block_size^2)
)
Out[ ]:
Mass Fraction -0.06 -0.05 -0.04 -0.03 -0.02 -0.01 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 0.11 -0.050 -0.045 -0.040 -0.035 -0.030 -0.025 -0.020 -0.015 -0.010 -0.005 0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035 0.040 0.045 0.050 0.055 0.060 0.065 0.070 0.075 0.080 0.085 0.090 0.095 0.100 0.105 -0.05 0.00 0.05 0.10 -0.050 -0.048 -0.046 -0.044 -0.042 -0.040 -0.038 -0.036 -0.034 -0.032 -0.030 -0.028 -0.026 -0.024 -0.022 -0.020 -0.018 -0.016 -0.014 -0.012 -0.010 -0.008 -0.006 -0.004 -0.002 0.000 0.002 0.004 0.006 0.008 0.010 0.012 0.014 0.016 0.018 0.020 0.022 0.024 0.026 0.028 0.030 0.032 0.034 0.036 0.038 0.040 0.042 0.044 0.046 0.048 0.050 0.052 0.054 0.056 0.058 0.060 0.062 0.064 0.066 0.068 0.070 0.072 0.074 0.076 0.078 0.080 0.082 0.084 0.086 0.088 0.090 0.092 0.094 0.096 0.098 0.100 0.102 Ca Ti Cr Mn Fe Ni Cu Zn Ag Ba C O Na Mg Al Si P S Cl K Element h,j,k,l,arrows,drag to pan i,o,+,-,scroll,shift-drag to zoom r,dbl-click to reset c for coordinates ? for help ? -1.80×10⁵ -1.60×10⁵ -1.40×10⁵ -1.20×10⁵ -1.00×10⁵ -8.00×10⁴ -6.00×10⁴ -4.00×10⁴ -2.00×10⁴ 0 2.00×10⁴ 4.00×10⁴ 6.00×10⁴ 8.00×10⁴ 1.00×10⁵ 1.20×10⁵ 1.40×10⁵ 1.60×10⁵ 1.80×10⁵ 2.00×10⁵ 2.20×10⁵ 2.40×10⁵ 2.60×10⁵ 2.80×10⁵ 3.00×10⁵ 3.20×10⁵ 3.40×10⁵ -1.60×10⁵ -1.50×10⁵ -1.40×10⁵ -1.30×10⁵ -1.20×10⁵ -1.10×10⁵ -1.00×10⁵ -9.00×10⁴ -8.00×10⁴ -7.00×10⁴ -6.00×10⁴ -5.00×10⁴ -4.00×10⁴ -3.00×10⁴ -2.00×10⁴ -1.00×10⁴ 0 1.00×10⁴ 2.00×10⁴ 3.00×10⁴ 4.00×10⁴ 5.00×10⁴ 6.00×10⁴ 7.00×10⁴ 8.00×10⁴ 9.00×10⁴ 1.00×10⁵ 1.10×10⁵ 1.20×10⁵ 1.30×10⁵ 1.40×10⁵ 1.50×10⁵ 1.60×10⁵ 1.70×10⁵ 1.80×10⁵ 1.90×10⁵ 2.00×10⁵ 2.10×10⁵ 2.20×10⁵ 2.30×10⁵ 2.40×10⁵ 2.50×10⁵ 2.60×10⁵ 2.70×10⁵ 2.80×10⁵ 2.90×10⁵ 3.00×10⁵ 3.10×10⁵ 3.20×10⁵ -2.0×10⁵ 0 2.0×10⁵ 4.0×10⁵ -1.600×10⁵ -1.550×10⁵ -1.500×10⁵ -1.450×10⁵ -1.400×10⁵ -1.350×10⁵ -1.300×10⁵ -1.250×10⁵ -1.200×10⁵ -1.150×10⁵ -1.100×10⁵ -1.050×10⁵ -1.000×10⁵ -9.500×10⁴ -9.000×10⁴ -8.500×10⁴ -8.000×10⁴ -7.500×10⁴ -7.000×10⁴ -6.500×10⁴ -6.000×10⁴ -5.500×10⁴ -5.000×10⁴ -4.500×10⁴ -4.000×10⁴ -3.500×10⁴ -3.000×10⁴ -2.500×10⁴ -2.000×10⁴ -1.500×10⁴ -1.000×10⁴ -5.000×10³ 0 5.000×10³ 1.000×10⁴ 1.500×10⁴ 2.000×10⁴ 2.500×10⁴ 3.000×10⁴ 3.500×10⁴ 4.000×10⁴ 4.500×10⁴ 5.000×10⁴ 5.500×10⁴ 6.000×10⁴ 6.500×10⁴ 7.000×10⁴ 7.500×10⁴ 8.000×10⁴ 8.500×10⁴ 9.000×10⁴ 9.500×10⁴ 1.000×10⁵ 1.050×10⁵ 1.100×10⁵ 1.150×10⁵ 1.200×10⁵ 1.250×10⁵ 1.300×10⁵ 1.350×10⁵ 1.400×10⁵ 1.450×10⁵ 1.500×10⁵ 1.550×10⁵ 1.600×10⁵ 1.650×10⁵ 1.700×10⁵ 1.750×10⁵ 1.800×10⁵ 1.850×10⁵ 1.900×10⁵ 1.950×10⁵ 2.000×10⁵ 2.050×10⁵ 2.100×10⁵ 2.150×10⁵ 2.200×10⁵ 2.250×10⁵ 2.300×10⁵ 2.350×10⁵ 2.400×10⁵ 2.450×10⁵ 2.500×10⁵ 2.550×10⁵ 2.600×10⁵ 2.650×10⁵ 2.700×10⁵ 2.750×10⁵ 2.800×10⁵ 2.850×10⁵ 2.900×10⁵ 2.950×10⁵ 3.000×10⁵ 3.050×10⁵ 3.100×10⁵ 3.150×10⁵ 3.200×10⁵ Counts

We can easily query the quantified data using broadcast syntax to create BitMatrix based on a boolean comparison. To determine which pixels have a mass-fraction of Si > 5% we create a mask. In this BitMatrix mask, 1 represents a true comparison and 0 a false comparison.

In [ ]:
mask_si = quant[n"Si"] .> 0.05
Out[ ]:
1024×1024 BitMatrix:
 0  0  0  1  0  1  1  1  1  0  0  0  0  …  0  0  0  0  0  0  0  0  1  1  0  0
 0  0  0  0  0  0  0  0  0  1  0  0  0     0  0  0  0  0  0  0  0  0  1  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  1  0  0
 0  0  0  0  0  0  0  0  0  0  0  1  0     0  0  0  0  0  1  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  1  0  0  0  0  1
 0  0  0  0  0  0  0  0  0  0  0  0  0  …  0  0  0  0  0  0  0  0  0  0  0  0
 0  1  1  0  0  0  0  0  0  0  0  0  0     1  0  0  0  0  0  0  0  0  0  0  0
 0  1  0  0  0  0  0  0  0  0  0  1  1     0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  1  1  0  0  0  0  0  1  1  1     0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  0  1  0  0  0  0  0  1  1  1     0  0  0  0  1  0  0  1  0  0  0  0
 0  0  0  0  0  1  0  0  0  0  0  1  1  …  0  0  0  0  0  0  0  1  1  1  0  0
 0  1  0  0  0  1  1  0  1  1  0  0  0     0  0  0  0  0  0  0  0  1  0  0  0
 0  0  0  0  0  1  1  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  0  0
 ⋮              ⋮              ⋮        ⋱           ⋮              ⋮        
 0  0  0  0  1  1  1  1  1  1  1  1  0     1  1  1  1  1  1  1  1  1  1  1  1
 0  0  0  1  0  1  1  1  1  1  1  1  1     1  1  1  1  1  1  1  1  0  0  1  1
 0  0  1  1  1  1  1  1  1  1  1  0  0     1  0  1  0  1  1  1  1  1  1  0  1
 1  1  1  0  0  1  1  1  1  0  1  1  0  …  1  0  0  1  0  0  1  1  1  1  0  1
 1  1  1  1  1  1  1  0  1  1  1  1  0     1  1  1  1  1  1  1  1  1  1  1  1
 1  1  1  1  1  1  0  0  0  0  1  0  1     1  1  1  1  0  1  1  1  1  1  0  0
 1  1  1  1  1  0  0  0  1  1  1  1  0     1  1  0  1  1  1  1  1  0  0  0  0
 0  1  1  0  0  0  0  0  0  0  0  0  0     1  1  1  1  1  1  1  1  1  1  0  0
 0  0  0  1  0  0  0  0  0  0  1  0  0  …  1  1  1  1  1  1  1  1  0  0  0  1
 0  1  1  0  0  0  0  0  0  0  0  0  0     0  1  0  1  1  0  1  1  1  1  1  1
 0  0  0  0  0  0  0  0  0  0  0  0  0     0  0  0  1  1  1  1  0  1  1  1  1
 0  0  0  1  1  1  0  0  0  0  0  0  0     0  0  0  0  0  1  0  0  1  1  1  1

Using similar syntax we can determine how many pixels have at least 50%, 10%, 5% and 1% of an element?

In [ ]:
DataFrame(
    Element=symbol.(els), 
    P50=map(el->count(quant[el] .> 0.5),els), 
    P10=map(el->count(quant[el] .> 0.1),els), 
    P5=map(el->count(quant[el] .> 0.05),els), 
    P1=map(el->count(quant[el] .> 0.01),els)
)
Out[ ]:
20×5 DataFrame
RowElementP50P10P5P1
StringInt64Int64Int64Int64
1C8670472085810199091045263
2O35574103499510406361045908
3Na0395515603161
4Mg096896727196
5Al0538039062742418
6Si5121297270774942104
7P02136532336
8S0128115253
9Cl00123493
10K04533768269937
11Ca011905723809920
12Ti01594694871
13Cr003121
14Mn184863934918021990386
15Fe33182310506314960971
16Ni21015543173
17Cu0611416359
18Zn27521485702
19Ag011217231
20Ba0212106467902

A mask can then be applied to extract and sum the requisite spectra within the hyperspectrum.

In [ ]:
plot( sum(hs, quant[n"Si"] .> 0.05), klms=els, xmax=10.0e3)
Out[ ]:
Energy (eV) -15000 -10000 -5000 0 5000 10000 15000 20000 25000 -10000 -9000 -8000 -7000 -6000 -5000 -4000 -3000 -2000 -1000 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 11000 12000 13000 14000 15000 16000 17000 18000 19000 20000 -10000 0 10000 20000 -10000 -9500 -9000 -8500 -8000 -7500 -7000 -6500 -6000 -5500 -5000 -4500 -4000 -3500 -3000 -2500 -2000 -1500 -1000 -500 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 5500 6000 6500 7000 7500 8000 8500 9000 9500 10000 10500 11000 11500 12000 12500 13000 13500 14000 14500 15000 15500 16000 16500 17000 17500 18000 18500 19000 19500 20000 Si Si K K K K Ti Ti Ti Ti Ag Ag Ag Ag Cu Cu Cu Mn Mn Mn Ag Ag Ag Ag Ag Cu Cu Al Al P P Mn Mn Mn Fe Fe Fe S S Ca Ca Ca Fe Fe Fe Mg Mg Cr Cr Cr Ni Ni Ni O O Cl Cl Zn Zn Zn K K Ca Ca Ca Ca C C Ti Ti Ti Na Na Ba Ba Ba Ba Ba Ba Ba Ba Ba Cr Cr Cr Ni Ni Zn Zn h,j,k,l,arrows,drag to pan i,o,+,-,scroll,shift-drag to zoom r,dbl-click to reset c for coordinates ? for help ? MaskedSum[Deep Sea Mn Nodule] Spectrum -30000000 -25000000 -20000000 -15000000 -10000000 -5000000 0 5000000 10000000 15000000 20000000 25000000 30000000 35000000 40000000 45000000 50000000 55000000 -26000000 -24000000 -22000000 -20000000 -18000000 -16000000 -14000000 -12000000 -10000000 -8000000 -6000000 -4000000 -2000000 0 2000000 4000000 6000000 8000000 10000000 12000000 14000000 16000000 18000000 20000000 22000000 24000000 26000000 28000000 30000000 32000000 34000000 36000000 38000000 40000000 42000000 44000000 46000000 48000000 50000000 -25000000 0 25000000 50000000 -25000000 -24000000 -23000000 -22000000 -21000000 -20000000 -19000000 -18000000 -17000000 -16000000 -15000000 -14000000 -13000000 -12000000 -11000000 -10000000 -9000000 -8000000 -7000000 -6000000 -5000000 -4000000 -3000000 -2000000 -1000000 0 1000000 2000000 3000000 4000000 5000000 6000000 7000000 8000000 9000000 10000000 11000000 12000000 13000000 14000000 15000000 16000000 17000000 18000000 19000000 20000000 21000000 22000000 23000000 24000000 25000000 26000000 27000000 28000000 29000000 30000000 31000000 32000000 33000000 34000000 35000000 36000000 37000000 38000000 39000000 40000000 41000000 42000000 43000000 44000000 45000000 46000000 47000000 48000000 49000000 50000000 Counts

Let's do something similar for Cu > 10%.

In [ ]:
plot( sum(hs, quant[n"Cu"] .> 0.1), klms=els, xmax=10.0e3)
Out[ ]:
Energy (eV) -15000 -10000 -5000 0 5000 10000 15000 20000 25000 -10000 -9000 -8000 -7000 -6000 -5000 -4000 -3000 -2000 -1000 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 11000 12000 13000 14000 15000 16000 17000 18000 19000 20000 -10000 0 10000 20000 -10000 -9500 -9000 -8500 -8000 -7500 -7000 -6500 -6000 -5500 -5000 -4500 -4000 -3500 -3000 -2500 -2000 -1500 -1000 -500 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 5500 6000 6500 7000 7500 8000 8500 9000 9500 10000 10500 11000 11500 12000 12500 13000 13500 14000 14500 15000 15500 16000 16500 17000 17500 18000 18500 19000 19500 20000 Si Si K K K K Ti Ti Ti Ti Ag Ag Ag Ag Cu Cu Cu Mn Mn Mn Ag Ag Ag Ag Ag Cu Cu Al Al P P Mn Mn Mn Fe Fe Fe S S Ca Ca Ca Fe Fe Fe Mg Mg Cr Cr Cr Ni Ni Ni O O Cl Cl Zn Zn Zn K K Ca Ca Ca Ca C C Ti Ti Ti Na Na Ba Ba Ba Ba Ba Ba Ba Ba Ba Cr Cr Cr Ni Ni Zn Zn h,j,k,l,arrows,drag to pan i,o,+,-,scroll,shift-drag to zoom r,dbl-click to reset c for coordinates ? for help ? MaskedSum[Deep Sea Mn Nodule] Spectrum -700 -600 -500 -400 -300 -200 -100 0 100 200 300 400 500 600 700 800 900 1000 1100 1200 1300 -600 -550 -500 -450 -400 -350 -300 -250 -200 -150 -100 -50 0 50 100 150 200 250 300 350 400 450 500 550 600 650 700 750 800 850 900 950 1000 1050 1100 1150 -1000 0 1000 2000 -580 -560 -540 -520 -500 -480 -460 -440 -420 -400 -380 -360 -340 -320 -300 -280 -260 -240 -220 -200 -180 -160 -140 -120 -100 -80 -60 -40 -20 0 20 40 60 80 100 120 140 160 180 200 220 240 260 280 300 320 340 360 380 400 420 440 460 480 500 520 540 560 580 600 620 640 660 680 700 720 740 760 780 800 820 840 860 880 900 920 940 960 980 1000 1020 1040 1060 1080 1100 1120 1140 Counts

Finally, to visualize the relative abundances and locations of the elements in the quantified data, we can construct a RGB image.

In [ ]:
using ImageCore
colorview(RGB, norm_quant[n"Mn"], norm_quant[n"Si"], norm_quant[n"O"])
Out[ ]:

Compare this with the ROI map from earlier.

In [ ]:
hs[n"Mn K-L3", n"Si K-L3", n"O K-L3"]
Out[ ]:

This is just a flavor of what you can do with spectra and hyperspectra using Julia, the Julia library infrastructure including the NeXL libraries. There are many machine learning, chemometric, image analysis and other libraries that can be readily applied.